Glassy Dislocation Dynamics in 2-D Colloidal Dimer Crystals 
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Although glassy relaxation is typically associated with disorder, here we report on a new type of 
glassy dynamics relating to dislocations within 2-D crystals of colloidal dimers. Previous studies have 
demonstrated that dislocation motion in dimer crystals is restricted by certain particle orientations. 
Here, we drag an optically trapped particle through such dimer crystals, creating dislocations. We 
find a two-stage relaxation response where initially dislocations glide until encountering particles that 
cage their motion. Subsequent relaxation occurs logarithmically slowly through a second process 
where dislocations hop between caged configurations. Finally, in simulations of sheared dimer 
crystals, the dislocation mean squared displacement displays a caging plateau typical of glassy 
dynamics. Together, these results reveal a novel glassy system within a colloidal crystal. 

PACS numbers: 82.70.Dd, 61.72.Ff, 64.70.pv 
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Dislocation mobility is central to both the materials 
properties and relaxation mechanisms of crystalline ma- 
terials [IHij]. Previous studies of dislocation motion in 
colloidal crystals composed of spherical particles have al- 
lowed a particle-scale view of defect formation and trans- 
port 0-0I- More recent experiments [8, 9] have begun to 
interrogate the role of particle anisotropy in determin- 
ing the rules of defect motion. In particular, studies of 
dislocation motion in crystals of colloidal dimer particles 
have uncovered novel restrictions on dislocation mobility. 
In these "degenerate crystals", the dimer lobes occupy 
triangular lattice sites while the particles are randomly 
oriented among the three crystalline directions, as shown 
in Fig. [TJ One consequence of the random orientations of 
the dimers is that dislocation glide is severely limited by 
certain particle arrangements in degenerate crystals Q. 

The present work utilizes local mechanical perturba- 
tion experiments to investigate the effects of this re- 
stricted dislocation motion on the relaxation mechanisms 
in degenerate crystals of colloidal dimers. Holographic 
optical tweezers are used to manipulate single lobe-sized 
spherical intruder particles within an otherwise pure de- 
generate crystal grain, deforming the crystal and intro- 
ducing defects. During the subsequent relaxation of the 
degenerate crystal, dislocations formed during the defor- 
mation leave the grain, either via annihilation with other 
dislocations or by moving to a grain boundary. Inter- 
estingly, we find that in large crystalline grains this dis- 
location relaxation occurs through a two-stage process 
reminiscent of slow relaxations in glassy systems, sug- 
gesting the novel concept that glassy phenomena may be 
present within certain kinds of colloidal crystals. 

We synthesize sterically stabilized silica dimer shell 
particles with spherical lobes of diameter 1.3 fim and 
lobe separation 1.4 /im as previously described [1, 0. 
For the spherical intruder particles, we include 1% vol- 



ume fraction of 1.3 /jm polystyrene spherical particles 
that are coated with silica and sterically stabilized us- 
ing polyvinylpyrrolidone (PVP) to produce surface chem- 
istry identical to the dimer particles. All particles are 
suspended in an aqueous solution of dimethylsulfoxide 
(DMSO) that index-matches the silica shells (with den- 
sity mismatch ~ 0.5 g/mL). The polystyrene cores of the 
spherical intruder particles remain index-mismatched, 
allowing for optical manipulation using laser tweezers. 
Thus, perturbations to the degenerate crystal are ap- 
plied only via the motion of the intruder particle. The 
suspension is pipetted into a sealed wedge-shaped glass 
cell as previously described Q and the particle area frac- 
tion (maintained at ~ 0.8) is controlled by tilting the cell 
so that particles sediment into the viewing region, which 
only accommodates a 2-D monolayer of particles. During 
the perturbation experiments, the crystal is imaged us- 
ing a confocal microscope integrated with the holographic 
optical tweezer system (Arryx, Inc) [loj . 

In preparation for a local perturbation experiment, one 
intruder particle is moved to the center of a degenerate 
crystal grain, and then the system is allowed to relax for 
at least 10 hours until the grain is defect-free. Fig. QJi-c 
shows a typical drag experiment for a small degenerate 
crystal grain consisting of N » 100 dimer particles. At 
time t = s (Fig. [Th.) the optical tweezers are turned on 
and used to drag the intruder particle by one lattice con- 
stant (LC) along a lattice direction, arriving at the new 
lattice position at t = 20 s. Then the optical tweezers 
are turned off and the grain is imaged until all the de- 
fects leave the grain or recombine (Fig. [TJ;) . Dislocations 
created by the drag deformation are identified from the 
microscope images by calculating the number of near- 
est neighbors for each lobe using Voronoi analysis. Each 
dislocation, consisting of paired five-fold and seven-fold 
coordinated lobes, is tracked over the experiment. 
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FIG. 1: Micrographs illustrating a local perturbation exper- 
iment. Five-fold and seven-fold coordinated defects created 
by the perturbation are marked with pink and blue dots, re- 
spectively, (a) A single spherical intruder particle in a small 
grain is dragged by 1 LC using optical tweezers, (b) Two 
dislocations marked by arrows glide to the grain edge and 
are absorbed by the grain boundary (solid closed curve), (c) 
The intruder particle remains stationary after the tweezers are 
turned off. The vacancy left behind has no topological charge 
and is stable, (d-f) In a large grain (only the central part is 
shown here) the dislocations marked by arrows are obstructed 
by the particles in red. When the tweezers are turned off, the 
blocked dislocations return along their glide paths, recombin- 
ing and causing the intruder particle to recoil. 



We observe different relaxation behavior depending on 
the degenerate crystal grain size. In small grains, the 
intruder particle remains stationary after being dragged 
to its new lattice position, even after the tweezers are 
turned off. In such cases, dislocations like those identi- 
fied by arrows in Fig. [Tb glide to the edge of the grain 
and are absorbed by the grain boundary, which we iden- 
tify using Voronoi analysis (details in the supplementary 
materials). In larger grains (Fig. [TJi-f ) , we observe that 
dislocations gliding toward the grain boundary are more 
likely to encounter glide-blocking particles whose orien- 
tations prevent further glide [8j. For example, the dislo- 



cations marked by arrows in Fig. must travel so far to 
reach the grain edge (beyond the field of view of Fig. QJ) 
that they encounter the blocking particles shown in red. 
Consequently, when the optical tweezers are turned off, 
these dislocations return toward the intruder particle, re- 
combining with the other defects produced by the defor- 
mation, and causing the intruder particle to recoil. 

To elucidate the effect of grain size on the recoil dis- 
tance of the intruder particle we conduct 19 independent 
deformation experiments on degenerate crystal grains of 
varying size. We find that a crossover from the station- 
ary response to the recoiling response occurs when dislo- 
cations produced during the deformation must traverse 
distances Z > 10 LC to reach a grain boundary (Fig. [5]). 
This result is consistent with the previously reported dis- 
tribution of unrestricted dislocation glide distances p(Z), 
which was found to decay exponentially (dashed gray 
curve in Fig. [2]) [8]. Together these results depict the 
probabilistic process, governed by the distribution p(Z), 
by which a deformation-induced dislocation may glide to 
a grain edge without meeting an obstacle. 

Although unrestricted dislocation glide over long dis- 
tances is improbable, multi-defect mechanisms such as 
dislocation reactions can allow dislocations to bypass 
glide-blocking obstacles and achieve long range motion 
[8J. However, since the energy required for such disloca- 
tion interactions is higher than that for simple dislocation 
glide [8]. larger mechanical deformations are required to 
access such processes. Furthermore, mechanisms involv- 
ing more than one defect are statistically less probable, 
and consequently may require longer waiting times. 
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FIG. 2: Recoil distance R as a function of glide distance Z to 
the nearest grain boundary. The error bars represent uncer- 
tainty in particle position. The distribution of unrestricted 
glide distances p(Z) (dashed gray curve) is reprinted from 



In order to probe such multi-defect relaxation re- 
sponses we conduct long time relaxation experiments 
where a perturbation is rapidly applied to a large degen- 
erate crystal grain. Before each experiment, one intruder 
particle is placed near the center of a large degenerate 
crystal grain containing more than 700 dimer particles. 
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At time t = s, the intruder particle is quickly dragged 
by 3 LC, arriving at its new lattice position at t = 5 s. 
The optical trap is turned off and the system is imaged 
for 10 hr, or until all the defects produced by the defor- 
mation leave the crystal grain. We record the number 
of defects by counting the lobes within the crystal grain 
that do not have six nearest neighbors. For example, a 
dislocation consisting of one five-fold and one seven-fold 
coordinated lobe is counted as two defects. 

To compare the relaxation response of degenerate crys- 
tals with a system where dislocations are known to glide 
without geometric restrictions, we repeat this experimen- 
tal procedure on crystals of spheres. Here we prepare 
crystals of PVP-stabilized 1 fim silica spheres suspended 
in an index-matching aqueous solution of DMSO, and 
we use 1 fim PVP-stabilized silica-coated polystyrene 
spheres as the intruder particles. 

The results from long time relaxation experiments 
on both degenerate crystals and crystals of spheres are 
shown in Fig. [3] In crystals of spheres, the decrease 
in the average number of defects Nd after the optical 
trap is turned off follows an exponential decay. The best 
fit exponential for this decay is Nd = 23e~*/ TS , where 
rs = 5 ± 0.5 s. In degenerate crystals, the initial decay 
of Nd is also characterized by an exponential decay, but 
a long tail is evident for later times. The best fit curve 
for this system combines both exponential and logarith- 
mic decay terms: Nd = 25e _t / r < 3 + 0.55 ln(l + ^f-), where 
Tp = 6 ± 1 s and r Q = 2 ± 1 x 10 5 s. These data depict 
strikingly different relaxation responses in the two sys- 
tems. In crystals of spheres, dislocations leave the grain 
via one fast mechanism with characteristic timescale ts- 
In degenerate crystals, a similar initial fast response with 
timescale rp is followed by a much slower process with a 
timescale r a that is 5 orders of magnitude larger than rg. 

Such two-stage relaxation is reminiscent of the slow re- 
laxation dynamics common to many glassy systems [lll - 
Il3j |. Studies of colloidal glasses have revealed a ca ging 
effect where particles are confined by their neighbors [1 3J 
and are transported through the glass in two stages: an 
initial fast diffusion until encountering the caging neigh- 
bors, and a much slower cage- hopping process caused by 
multi-particle rearrangements that shift the cage struc- 
ture. This two-stage glassy relaxation is characterized by 
a long tail similar to that observed in defect relaxation 
for degenerate crystals (Fig. [3]), suggesting that defect 
dynamics in dimer crystals is glassy. 

Similarly to a particle in a glass, a dislocation in a 
degenerate crystal is caged. Initially, it moves via un- 
restricted glide at a rate characterized by Tg, until it 
reaches particles with glide-blocking orientations. Here it 
remains caged until a second process characterized by the 
time scale r Q can allow it to hop to a new region of unre- 
stricted glide. This two stage relaxation process also ex- 
plains the recoil measurements in Fig. [5] For small grains 
where dislocations only need to travel a short distance to 
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FIG. 3: Average number of defects versus time for large 
(> 700 particles) degenerate crystals and crystals of spheres. 
Error bars represent the standard error of the mean. The solid 
line is the best fit exponential for the decay in the sphere data 
(R 2 — 0.99); the dashed line is the best fit sum of an exponen- 
tial and a logarithmic decay for the dimer data (R 2 = 0.99). 



reach the grain boundary, relaxation can occur through 
unrestricted glide alone. For larger grains, both stages 
of relaxation would be needed for the dislocations to exit 
the grain. However, this would require sufficiently long 
time scales for cage hopping processes to occur. Particles 
in the recoil experiments are only dragged 1 LC with a 
drag duration of t 3 < 20 s <C r Q , which does not create 
enough defects to trap the system and allow enough time 
for multidefect relaxation processes, thus leading to the 
observed recoil in large grains. 

Glassy relaxation can be directly observed in colloidal 
glasses by measuring the particle mean square displace- 
ment (MSD) versus time. Visualizing dislocation trajec- 
tories in local deformation experiments is complicated by 
the fact that all the defects originate near the intruder 
particle, rather than uniformly throughout the system 
as would occur during macroscopic deformations. To ex- 
plore such macroscopic deformations we run Non Equilib- 
rium Molecular Dynamics (NEMD) simulations of dimer 
and sphere crystals under uniform shear. This approach 
probes defect dynamics in a non-equilibrium "steady 
state" rather than a perturbation. Furthermore, such 
simulations allow for comparison of dimer and sphere 
crystals in ideal conditions without grain boundaries or 
particle polydispersity. The simulations are conducted 
in a 2-D canonical ensemble wherein hard spheres (or 
dimer lobes) interact via the Weeks-Chandler- Andersen 
potential |l4|. In reduced Lennard- Jones units, the sys- 
tem has size 50 x 50, sphere or lobe diameter d = 1, 
dimer lobe separation 1.07 maintained using holonomic 
constraints, and temperature of 1 maintained via a con- 
figurational thermostat [HI]. Area fractions of 0.808 for 
dimers and 0.792 for spheres are chosen to have both 
systems at identical pressure, temperature, and chemical 
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simulation time steps (St) 

FIG. 4: (Color online) Defect trajectories in simulated crys- 
tals under shear (large arrows), (a) Dislocations in a crystal 
of spheres glide along straight paths in response to the applied 
shear, (b) Dislocations in a degenerate crystal follow crooked 
trajectories as they hop between caged configurations, (c) 
MSD of individual dislocations in crystals of spheres and de- 
generate crystals. Error bars are smaller than the markers. 
The solid curve is the best fit for MSD = ( Ar 2 ) = u 2 t 2 + 2Dt. 



ical of particles in glassy systems. 

In conclusion, we have uncovered a novel glassy system 
where the constituent particles are assembled into an or- 
dered crystalline structure, but the dislocations within 
these crystals are caged and demonstrate slow, two-stage 
glassy relaxation. Such crystals also have a natural mech- 
anism for encoding memory effects commonly observed 
in glasses, since moving dislocations reorient the dimcr 
particles, encoding constraints for subsequent defect mo- 
tion. While the transition into the glassy state typ ically 
occurs through an increase in particle density [18j . here 
the glassy dislocation dynamics arises from constraints 
in an already dense crystal. Future studies of crystals of 
spheres doped with increasing concentrations of dimers 
will probe for this transition and determine whether it 
can be thought of more gcncrically as an additional route 
to the jammed state. Finally, unlike typical jammed sys- 
tems, increasing the dislocation density can reduce the 
amount of caging, suggesting that degenerate crystals 
may be self-healing materials. 

We thank R. Ganapathy, X. Cheng, and the A. Liu 
group for helpful discussion. This work was supported 
by the DOE, Basic Energy Sciences, Grant #ER46517. 



potential. A homogeneous shear field is applied to ini- 
tially defect-free crystals using SLLOD equations of mo- 
tion [161 ] and Lees-Edwards periodic boundary conditions 
[l7| . The equations of motion are integrated using the 
Runge Kutta 4th order method with time step St = 0.01. 
Uniform shear is applied at a strain rate 7 = 5xl0~ 5 per 
St, a value small enough to drive dislocation motion at a 
speed independent of 7. Dislocations are monitored af- 
ter the number of defects has reached a steady value (full 
simulation details provided in the supplementary mate- 
rials section, below). 

These simulations confirm strikingly different defect 
transport mechanisms in the two systems. Fig. 0^-b 
show plots of defect positions, where color corresponds 
to time. In crystals of spheres, dislocation pairs form 
and glide apart along straight paths nearly aligned with 
the shear direction (Fig. In stark contrast, de- 

fects in degenerate crystals follow crooked trajectories 
(Fig. [3J5) and do not glide throughout the crystal, but 
rather are trapped within local cages. This caging is 
evident in the dislocation MSD, shown for both degen- 
erate crystals and crystals of spheres in Fig. To cal- 
culate the dislocation MSD versus time, individual dis- 
locations that can be resolved from other nearby de- 
fects are tracked between time steps. In crystals of 
spheres, dislocation motion is consistent with a biased 
ID random walk (solid curve, Fig. 0J;) with a diffusion 
constant D = 2.3 ± 0.2 x 10~ 2 d 2 /St, and drift speed 
u = 2.23 ±0.05 x 10" 2 d/St. The dislocation MSD in de- 
generate crystals displays the characteristic plateau typ- 
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GLASSY DISLOCATION DYNAMICS IN 2-D 
COLLOIDAL DIMER CRYSTALS: 
SUPPLEMENTARY MATERIALS 

MOLECULAR DYNAMICS SIMULATIONS 

We have performed 2-D Number Volume Temperature 
(NVT) Non Equilibrium Molecular Dynamics (NEMD) 
simulations of both pure sphere and pure dimer crystals. 
The separation between lobes in the dimer particles is 
set to be 1.07 times the diameter of the spherical lobes, 
in approximation of the experimental system. The hard 
spheres (or dimer lobes) interact with each other via the 
Weeks-Chandler- Andersen potential given by 

u(r) = f 4(r~ 12 - r- 6 ) + 1 r < r m = 2* 
\ r>r m 

where r is the interparticle distance and r m is the mini- 
mum of the potential. In this and all the following equa- 
tions, all the properties are expressed in terms of reduced 
Lennard-Jones units. 

The crystals are initially defect-free, with every sphere 
(or dimer lobe) located exactly on a triangular lattice 
site. Dimer crystals are in a degenerate crystalline state 
such that the dimers are oriented randomly along any of 
the three lattice directions [19j] . The simulations of both 
dimer and sphere crystals are performed at equivalent 
area fractions, as defined below, at reduced temperatures 
T* = 1.0 (sufficiently below their melting temperatures) 
with a time step St = 10~ 2 . To preserve the hetero- 
geneities in the system we use a novel configurational 
thermostat which does not constrain the flow profile 

Area fraction 

In order to fairly compare crystals of dimers with crys- 
tals of spheres at equivalent area fractions, we first per- 
form direct interfacial NVT simulations of compound sys- 
tems of sphere and dimer lattices. On equilibration, once 
both lattices share the same chemical potential, pressure 
and temperature, these simulations provide estimates of 
the coexistence area fractions for the two species. Then, 
for our NEMD simulations, we have used one set of these 
equivalent area fractions: 0.808 for dimers and 0.792 for 
spheres, both well above the melting densities. 

Shear rate 

Both the crystals of pure spheres and the crystals of 
pure dimers are sheared at a shear rate 7 = 5 x 10~ 5 . 
In order to choose the value of 7, several factors were 
considered. First, in order to study dislocation motion 
in the crystalline regime, the shear rate must not be high 
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FIG. 5: Dislocation mean squared displacement in crystals of 
spheres under uniform shear at three different shear rates 7. 
Error bars represent the standard error of the mean. 

enough to destroy the crystalline order. When shear is 
applied to the initially defect-free system, some defects 
are created, thereby reducing the crystalline order. How- 
ever, if the shear rate is chosen to be low enough, then 
after some characteristic 'stabilization time' - about 10 6 
time steps in our simulations - the number of defects 
in the system fluctuates around a steady value and the 
crystalline order does not decrease further. 

A second important consideration is the effect of shear 
rate on dislocation glide speed in crystals of spheres. 
To ensure that we are studying dislocation motion in 
a 'quasi-static' regime where dislocation motion is inde- 
pendent of the applied shear rate, we have also tracked 
dislocation motion in crystals sheared at two additional 
shear rates: 7 — 25 x 10~ 5 and 7 = 1 x 10~ 5 . We com- 
pare the dislocation mean squared displacement (MSD) 
for the 3 shear rates and find that while dislocation mo- 
tion is slightly faster for 7 = 25 x 10~ 5 , the lower two 
shear rates are indistinguishable within the error bars 
(Fig. [5]). This indicates that at 7 = 5 x 10~ 5 , we are in 
the regime where dislocation motion is independent of 7. 

This shear rate independence can be understood by 
considering dislocation motion in crystals of spheres, 
where we find that long periods of dislocation inactiv- 
ity are punctuated by brief glide events. In such glide 
events, two dislocations form and glide apart along the 
lattice axis aligned with the shear direction. Between the 
glide events, shear stress caused by the applied shear rate 
7 slowly builds up until enough stress is present to sep- 
arate a pair of dislocations. Once the dislocations start 
to separate, they quickly glide apart and annihilate after 
traversing the periodic boundary conditions. Changing 
the shear rate 7 only affects the average waiting time be- 
tween glide events, not the speed of dislocation motion 
during the events. To quantify how small the applied 
shear rate is, we compare the time r~ t for the applied shear 
to cause a shear displacement of 1 LC to the time Td for 
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FIG. 6: The grain boundary defects are highlighted around 
the small grain also shown in Fig. lc of the manuscript. Pink 
and blue dots represent lobes having a number of nearest 
neighbors equal to 5 or 7, respectively. The dotted circle 
marks the spherical intruder particle. 



a dislocation pair to separate across the simulation box, 
causing a shear displacement of 1 LC. For j — 5 x 10~ 5 , 
r-y = 4 x 10 4 , while = 990 ± 16 (mean ± standard error 
of the mean). Thus, is about 40 times larger than r^. 



Calculating dislocation mean squared displacement 

To calculate the dislocation MSD, we have tracked the 
positions of individual dislocations consisting of paired 
5-fold and 7- fold coordinated lobes (or spheres). In both 
degenerate crystals of dimers and crystals of spheres, our 
analysis is restricted to dislocations that can be distinctly 
resolved from any other nearby defects. In dimer crystals, 
this requirement limits our analysis to dislocation motion 
before or after - but not during - a dislocation reaction. 

In order to correct for affine displacements imposed 
by the applied external shear 7, relative dislocation x- 
displacements Ai(t) were shifted by the accumulated 
strain displacement Sx(t) corresponding to the relative 
^-displacement Ay(r). The shear-corrected dislocation 
displacement Ax c (t) after a time interval r is given by: 

Ax c (t) = Ax(t)-Sx(t) 

= Ax(t) — 7rAy(r) 

Typical values of Sx(t) for time intervals r measured in 
calculating the dislocation MSD were approximately 3 
orders of magnitude smaller than Ax c (t). 



EXPERIMENTAL TECHNIQUES 



Computing grain boundaries 



To locate the grain boundaries in the experimental de- 
generate crystals, Voronoi analysis is used to compute 
the number of nearest neighbors for each lobe. Crys- 
talline grains containing many ordered lobes with 6 near- 
est neighbors each are surrounded by a string-like collec- 
tion of defect lobes having a number of nearest neighbors 
not equal to 6. The grain boundary defects for the small 
grain shown in Fig. lc of the manuscript are explicitly 
highlighted in Fig. [6] 



[1 
[2. 
[3 
[4 

[5 

[6 

[7 

[8] 
[9] 

[10; 
[11 
[12 

[13 
[14 

[15; 
[16 

[17 
[18 
[19 



E. Orowan, Z. Phys. 89, 605 (1934). 
M. Polanyi, Z. Phys. 89, 660 (1934). 

G. I. Taylor, Proc. R. Soc. A 145, 362 (1934). 

P. Schall, I. Cohen, D. Weitz, and F. Spaepen, Science 

305, 1944 (2004). 

A. Libal, C. Reichhardt, and C. J. O. Reichhardt, Phys. 
Rev. E 75, 011403 (2007). 

K. Zahn, R. Lenke, and G. Maret, Phys. Rev. Lett. 82, 
2721 (1999). 

V. de Villeneuve, R. Dullens, D. Aarts, E. Groeneveld, 
J. Scherff, W. Kegel, and H. Lekkerkerker, Science 309, 
1231 (2005). 

S. J. Gerbode, S. H. Lee, C. M. Liddell, and I. Cohen, 

Phys. Rev. Lett. 101, 058302 (2008). 

S. H. Lee, S. J. Gerbode, B. S. John, A. K. Wolfgang, 

F. A. Escobedo, I. Cohen, and C. M. Liddell, J. Mater. 
Chem. 18, 4912 (2008). 

E. R. Dufresne and D. G. Grier, Rev. Sci. Instr. 69, 1974 
(1998). 

W. Gotze and L. Sjogren, J. Phys.: Condens. Matter 1, 
4203 (1989). 

A. Coniglio, A. Fierro, H. J. Herrmann, and 
M. Nicodemi, eds., Unifying Concepts in Granular Media 
and Glasses (Elsevier, 2004). 

E. R. Weeks and D. A. Weitz, Phys. Rev. Lett. 89, 
095704 (2002). 

J. D. Weeks, D. Chandler, and H. C. Andersen, J. Chem. 
Phys. 54, 5237 (1971). 

C. Braga and K. P. Travis, J. Chem. Phys. 123, 134101 
(2005). 

D. J. Evans and G. P. Morriss, Statistical Mechanics of 
Nonequilibrium Liquids (Academic Press, 1990). 

A. W. Lees and S. F. Edwards, J. Phys. C 5, 1921 (1972). 
A. J. Liu and S. R. Nagel, Nature 396, 21 (1998). 
K. W. Wojciechowski, D. Frenkel, and A. C. Branka, 
Phys. Rev. Lett. 66, 3168 (1991). 



